7-2 WHO 8-point Outcome Scale

Working analyses of the WHO 8-point outcome scale outcome.

Authors
Affiliations
James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Published

July 13, 2022

Preamble

Load packages
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(ggdist)
library(patchwork)
library(lubridate)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))
Code
logit <- function(x) log(x) - log(1 - x)
expit <- function(x) 1 / (1 + exp(-x))
ordered_logit <- function(x) {
  c(
    1 - expit(-x[1]),
    expit(-x[1:(length(x)-1)]) - expit(-x[2:(length(x))]),
    expit(-tail(x, 1))
  )
}
Prepare dataset
devtools::load_all()
all_dat <- read_all_no_daily()

acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()
acs_itt_nona_dat <- acs_itt_dat %>%
  filter(!is.na(D28_who))

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(D28_who))

# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(D28_who))

# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(D28_who))
Load models
ordmod0 <- cmdstan_model(
  "stan/ordinal/logistic_cumulative.stan") # No epoch or site
ordmod_epoch <- cmdstan_model(
  "stan/ordinal/logistic_cumulative_epoch.stan") # Epoch only
ordmod <- cmdstan_model(
  "stan/ordinal/logistic_cumulative_epoch_site.stan") # Full model
logistic <- cmdstan_model(file.path("stan", "binary", "logistic_site_epoch.stan"))
Helper functions
make_summary_table <- function(dat, format = "html") {
  tdat <- dat %>%
  group_by(CAssignment = factor(CAssignment, 
                                levels = c("C1", "C2", "C3", "C4"),
                                labels = intervention_labels2()$CAssignment[-1])) %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(D28_who)),
    Deaths = sprintf("%i (%.0f%%)", 
                     sum(D28_who == 8, na.rm = TRUE), 
                     100 * mean(D28_who == 8, na.rm = TRUE)),
    `Hospitalised` = sprintf("%i (%.0f%%)", 
                     sum(D28_who > 2 & D28_who < 8, na.rm = TRUE), 
                     100 * mean(D28_who > 2 & D28_who < 8, na.rm = TRUE)),
    `WHO, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(D28_who, na.rm = T), 
      quantile(D28_who, 0.25, na.rm = TRUE), 
      quantile(D28_who, 0.75, na.rm = TRUE))
  ) %>%
  bind_rows(
    dat %>%
  group_by(CAssignment = "Overall") %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(D28_who)),
    Deaths = sprintf("%i (%.0f%%)", 
                     sum(D28_who == 8, na.rm = TRUE), 
                     100 * mean(D28_who == 8, na.rm = TRUE)),
    `Hospitalised` = sprintf("%i (%.0f%%)", 
                     sum(D28_who > 2 & D28_who < 8, na.rm = TRUE), 
                     100 * mean(D28_who > 2 & D28_who < 8, na.rm = TRUE)),
    `WHO, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(D28_who, na.rm = T), 
      quantile(D28_who, 0.25, na.rm = TRUE), 
      quantile(D28_who, 0.75, na.rm = TRUE))
  )
  ) %>%
  rename(`Anticoagulation\nintervention` = CAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}

make_primary_model_data <- function(
    dat, 
    vars = NULL,
    beta_sd_var = NULL,
    ctr = contr.orthonorm,
    outcome = "D28_who",
    p_mult = 2) {
  
  X <- make_X_design(dat, vars, ctr)
  attX <- attr(X, "contrasts")
  X <- X[, -1]
  attr(X, "contrasts") <- attX
  nXtrt <- sum(grepl("rand", colnames(X)))
  
  beta_sd <- c(rep(1, nXtrt), beta_sd_var)
  
  epoch <- dat$epoch
  M_epoch <- max(dat$epoch)
  region <- dat$ctry_num
  M_region <- max(region)
  site <- dat$site_num
  M_site <- max(site)
  region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)

  y_raw <- ordered(dat[[outcome]])
  y_mod <- as.integer(y_raw)
  N <- dim(X)[1]
  K <- dim(X)[2]
  unique_y <- length(unique(y_mod))
  list(N = N, K = K, J = unique_y, X = X, 
       y = y_mod, y_raw = y_raw,
       M_region = M_region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       p_par = p_mult * rep(1 / unique_y, unique_y),
       beta_sd = beta_sd)  
}

Descriptive

Anticoagulation

Summary of WHO outcome by arm
save_tex_table(
  make_summary_table(acs_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-2-anticoagulation-summary"))
make_summary_table(acs_itt_dat)
Anticoagulation intervention Patients Known Deaths Hospitalised WHO, Median (Q1, Q3)
Low-dose 610 596 19 (3%) 7 (1%) 1 (1, 2)
Intermediate-dose 613 603 15 (2%) 5 (1%) 1 (1, 2)
Low-dose with aspirin 283 281 10 (4%) 5 (2%) 1 (1, 2)
Therapeutic-dose 50 50 6 (12%) 1 (2%) 1 (1, 2)
Overall 1556 1530 50 (3%) 18 (1%) 1 (1, 2)
Table 1: Summary of WHO scale at day 28 by treatment group.
Code
p <- acs_itt_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", "\n")), 
    D28_who = ordered(as.integer(D28_who), levels = 1:8)
  ) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(CAssignment, p)) +
  geom_col(aes(fill = D28_who)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "WHO scale", y = "Cumulative proportion", x = "Anticoagulation intervention")
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-2-descriptive.pdf"), p, height = 3, width = 6)
p

Code
acs_itt_nona_dat %>%
  dplyr::count(CAssignment, D28_who = factor(as.integer(D28_who))) %>%
  complete(CAssignment, D28_who, fill = list(n = 0)) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(D28_who, p)) +
  facet_wrap( ~ CAssignment) +
  geom_point() +
  geom_segment(aes(xend = D28_who, y = 0, yend = p))

Age

Code
pdat <- all_dat %>%
  filter_acs_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    agegrp = cut(AgeAtEntry, c(18, seq(25, 75, 5), 100), include.lowest = T, right = F),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(agegrp) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(agegrp) %>%
  summarise(n = sum(n))

p1 <- ggplot(pdat2, aes(agegrp, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry") +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45))

p2 <- ggplot(pdat, aes(agegrp, p)) +
  geom_col(aes(fill = who)) +
  labs(x = "Age", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-age.pdf"), p, height = 2.5, width = 6)
p

Country

Code
pdat <- all_dat %>%
  filter_acs_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(Country) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
pdat2 <- pdat %>%
  group_by(Country) %>%
  summarise(n = sum(n))

p1 <- ggplot(pdat2, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")

p2 <- ggplot(pdat, aes(Country, p)) +
  geom_col(aes(fill = who)) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-country.pdf"), p, height = 2.5, width = 6)
p

Site

Code
pdat <- all_dat %>%
  filter_acs_itt() %>%
  filter(!is.na(D28_who)) %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    Site = fct_infreq(Location),
    who = ordered(D28_who, levels = 1:8)) %>%
  complete(who = ordered(1:8), nesting(Country, Site), fill = list(n = 0)) %>%
  group_by(Country, Site) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup() %>%
  mutate(
    Country = droplevels(Country),
    Site = droplevels(Site)
  )
pdat2 <- pdat %>%
  group_by(Country, Site) %>%
  summarise(n = sum(n)) %>%
  ungroup()
p1 <- ggplot(pdat2, aes(Site, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(pdat, aes(Site, p)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col(aes(fill = who)) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F, ncol = 1)) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 / p2 +
  plot_layout(guides = 'collect')
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-country-site.pdf"), p, height = 4, width = 6.25)
p

Calendar Time

Code
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    yr = year(RandDate), mth = month(RandDate),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(yr, mth) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p1 <- pdat %>%
  group_by(yr, mth) %>%
  summarise(n = sum(n)) %>%
  ggplot(., aes(mth, n))  +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p2 <- ggplot(pdat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
  geom_col(aes(fill = who)) +
  labs(x = "Calendar date (month of year)", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines")) +
  scale_x_continuous(breaks = 1:12)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-2-calendar-time.pdf"), p, height = 2, width = 6)
p

Sample Cumulative Logits

Proportional odds looks reasonable for all logits except 1.

Plot sample cumulative logits
trt_counts <- acs_itt_nona_dat %>%
  dplyr::count(CAssignment, D28_who) %>%
  complete(CAssignment, D28_who, fill = list(n = 0)) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(CAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(D28_who) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
  filter(D28_who != 28)

ggplot(trt_logit, aes(CAssignment, rel_clogit)) +
  facet_wrap( ~ D28_who) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Plot sample cumulative logits
ggplot(trt_logit, aes(D28_who, rel_clogit)) +
  facet_wrap( ~ CAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Modelling

Primary Model

Fit primary model
fit_primary_model <- function(dat, ctr = contr.orthonorm, save = FALSE) {
  mdat <- make_primary_model_data(
    dat, c("inelgc3", "agegte60", "ctry"), 
    c(10, 2.5, 1, 1), 
    ctr)
  snk <- capture.output(
    mfit <- ordmod$sample(
      data = mdat,
      seed = 712508,
      chains = 8,
      parallel_chains = 8,
      iter_warmup = 1000,
      iter_sampling = 2500,
      refresh = 0,
      adapt_delta = 0.98
  ))
  if (save) {
    mfit$save_object(file.path("outputs", "models", "secondary", "7-2.rds"))  
  }
  mdrws <- as_draws_rvars(mfit$draws())
  # Add names
  epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)
  site_map <- dat %>% dplyr::count(site_num, site)
  names(mdrws$beta) <- colnames(mdat$X)
  names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  names(mdrws$gamma_site) <- site_map$site
  # Convert to treatment log odds ratios
  mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
  mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
  mdrws$trtA <- mdrws$Acon[-1] - mdrws$Acon[1]
  mdrws$trtC <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  mdrws$compare <- c(
    "Intermediate vs low" = exp(mdrws$trtC[1]),
    "Low with aspirin vs low" = exp(mdrws$trtC[2]),
    "Therapeutic vs low" = exp(mdrws$trtC[3]),
    "Intermediate vs low with aspirin" = exp(mdrws$trtC[1] - mdrws$trtC[2]),
    "Intermediate vs therapeutic" = exp(mdrws$trtC[1] - mdrws$trtC[3]),
    "Low with aspirin vs therapeutic" = exp(mdrws$trtC[2] - mdrws$trtC[3])
  )
  mdrws$OR <- c(
    setNames(exp(mdrws$trtC), 
             c("Intermediate", "Low with aspirin", "Therapeutic")),
    "Ineligible aspirin" = exp(mdrws$beta[which(grepl("inelg", colnames(mdat$X)))]),
    "Age \u2265 60" = exp(mdrws$beta[which(grepl("age", colnames(mdat$X)))]),
    setNames(exp(mdrws$beta[which(grepl("ctry", colnames(mdat$X)))]), 
             c("Australia/New Zealand", "Nepal"))
  )
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}
res <- fit_primary_model(acs_itt_nona_dat, save = TRUE)
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/secondary/7-2-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.81 (0.61, 1.07) 0.82 (0.12) 0.93
Low with aspirin 0.76 (0.53, 1.09) 0.77 (0.14) 0.93
Therapeutic 1.53 (0.77, 3.06) 1.63 (0.60) 0.12
Ineligible aspirin 1.72 (0.78, 3.77) 1.86 (0.77) 0.09
Age ≥ 60 2.47 (1.90, 3.27) 2.50 (0.35) 0.00
Australia/New Zealand 1.37 (0.42, 4.25) 1.62 (1.02) 0.30
Nepal 0.73 (0.23, 2.46) 0.88 (0.63) 0.71
Posterior summaries for model parameters (fixed-effects).
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch),
  exp(res$drws$gamma_site),
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-2-primary-model-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Code
p <- plot_or_densities(res$drws$compare)
p

Code
res$fit$summary()
# A tibble: 3,173 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -1.06e+3 -1.06e+3 6.87e+0 6.78e+0 -1.08e+3 -1.05e+3  1.00    4798.
 2 p0[1]       8.08e-1  8.21e-1 7.83e-2 7.17e-2  6.61e-1  9.11e-1  1.00    4178.
 3 p0[2]       1.71e-1  1.60e-1 6.72e-2 6.32e-2  8.00e-2  2.97e-1  1.00    4205.
 4 p0[3]       1.02e-3  7.95e-4 8.47e-4 5.71e-4  2.01e-4  2.56e-3  1.00    7625.
 5 p0[4]       1.63e-3  1.33e-3 1.17e-3 8.34e-4  4.38e-4  3.84e-3  1.00    6787.
 6 p0[5]       1.32e-3  1.05e-3 1.01e-3 7.02e-4  3.12e-4  3.17e-3  1.00    6937.
 7 p0[6]       1.30e-3  1.04e-3 9.73e-4 6.93e-4  3.09e-4  3.16e-3  1.00    7225.
 8 p0[7]       6.90e-4  5.08e-4 6.36e-4 4.09e-4  9.84e-5  1.91e-3  1.00    9201.
 9 p0[8]       1.50e-2  1.31e-2 8.53e-3 6.48e-3  5.64e-3  3.10e-2  1.00    4369.
10 beta_raw[…  5.04e-1  5.05e-1 2.76e-1 2.76e-1  5.16e-2  9.60e-1  1.00   14961.
# … with 3,163 more rows, and 1 more variable: ess_tail <dbl>
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 1 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7457517 0.7515987 0.8222353 0.7093861 0.7255039 0.7288769 0.7351114
[8] 0.7027627
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["alpha"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["gamma_epoch"])

Posterior Predictive

Code
y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc1 <- function(grp) {
  ppc_dat  %>%
  group_by(grp = !!grp) %>%
  summarise(
    y_1 = mean(D28_who == 1),
    ypp_1 = rvar_mean(y_ppc == 1),
    y_2 = mean(D28_who <= 2),
    ypp_2 = rvar_mean(y_ppc <= 2),
    y_3 = mean(D28_who <= 3),
    ypp_3 = rvar_mean(y_ppc <= 3),
    y_4 = mean(D28_who <= 4),
    ypp_4 = rvar_mean(y_ppc <= 4),
    y_5 = mean(D28_who <= 5),
    ypp_5 = rvar_mean(y_ppc <= 5),
    y_6 = mean(D28_who <= 6),
    ypp_6 = rvar_mean(y_ppc <= 6),
    y_7 = mean(D28_who <= 7),
    ypp_7 = rvar_mean(y_ppc <= 7),
    y_8 = mean(D28_who <= 8),
    ypp_8 = rvar_mean(y_ppc <= 8)
  ) %>%
  pivot_longer(y_1:ypp_8, names_to = c("response", "who"), names_sep = "_", values_to = "posterior") %>%
  mutate(who = as.numeric(who))
}

grp_ppc2 <- function(grp) {
  ppc_dat  %>%
  group_by(grp = !!grp) %>%
  summarise(
    y_1 = mean(D28_who == 1),
    ypp_1 = rvar_mean(y_ppc == 1),
    y_2 = mean(D28_who == 2),
    ypp_2 = rvar_mean(y_ppc == 2),
    y_3 = mean(D28_who == 3),
    ypp_3 = rvar_mean(y_ppc == 3),
    y_4 = mean(D28_who == 4),
    ypp_4 = rvar_mean(y_ppc == 4),
    y_5 = mean(D28_who == 5),
    ypp_5 = rvar_mean(y_ppc == 5),
    y_6 = mean(D28_who == 6),
    ypp_6 = rvar_mean(y_ppc == 6),
    y_7 = mean(D28_who == 7),
    ypp_7 = rvar_mean(y_ppc == 7),
    y_8 = mean(D28_who == 8),
    ypp_8 = rvar_mean(y_ppc == 8)
  ) %>%
  pivot_longer(y_1:ypp_8, names_to = c("response", "who"), 
               names_sep = "_", values_to = "posterior") %>%
  mutate(who = as.numeric(who))
}

plot_grp_ppc <- function(dat) {
  ggplot(dat %>% filter(response == "ypp"), aes(x = who)) +
    facet_wrap( ~ grp, nrow = 1) +
    stat_slabinterval(aes(ydist = posterior))  +
    geom_point(data = dat %>% filter(response == "y"), 
               aes(x = who, y = mean(posterior)),
               colour = "red",
               shape = 23) +
    labs(x = "WHO outcome", 
         y = "Probability")
}

plot_grp_ppc <- function(dat, lab = "", xlab = "Probability") {
  ggplot(dat %>% filter(response == "ypp"), aes(y = who)) +
    facet_wrap( ~ grp, nrow = 1) +
    stat_slabinterval(aes(xdist = posterior), fatten_point = 1)  +
    geom_point(data = dat %>% filter(response == "y"), 
               aes(y = who, x = mean(posterior)),
               colour = "red",
               shape = 23) +
    scale_x_continuous(xlab, breaks = c(0, 0.5),
                       sec.axis = sec_axis(~ . , name = lab, breaks = NULL, labels = NULL)) +
    scale_y_continuous("WHO\noutcome", breaks = c(1,3,5,7)) +
    theme(strip.text = element_text(size = rel(0.7)),
          axis.title.x = element_text(size = rel(0.7)),
          axis.text.x = element_text(size = rel(0.65)),
          axis.title.y = element_text(size = rel(0.75)),
          axis.title.x.bottom = element_blank())
}

pp_epoch <- grp_ppc2(sym("epoch")) %>% 
  mutate(grp = fct_inorder(factor(grp)))
pp_C <- grp_ppc2(sym("CAssignment"))
pp_ctry <- grp_ppc2(sym("country"))
pp_site <- grp_ppc2(sym("site")) %>%
  left_join(ppc_dat %>% dplyr::count(site, country), by = c("grp" = "site"))

p1 <- plot_grp_ppc(pp_C, "Anticoagulation", "")
p2 <- plot_grp_ppc(pp_ctry, "Country", "") 
p3 <- plot_grp_ppc(pp_epoch, "Epoch", "")
p4 <- plot_grp_ppc(pp_site %>% filter(country == "IN"), "Sites India", "")
p5 <- plot_grp_ppc(pp_site %>% filter(country == "AU"), "Sites Australia", "")
p6 <- plot_grp_ppc(pp_site %>% filter(country == "NP"), "Sites Nepal", "")
p7 <- plot_grp_ppc(pp_site %>% filter(country == "NP"), "Sites New Zealand", "")
p <- (p1 | p2) / p3 / p4 / p5 / (p6 | p7)
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-2-primary-model-acs-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.75)
p

Sensitivity: Concurrent controls

Intermediate dose

  • Set: ACS-ITT-intermediate
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc2_dat, "latex"), 
  file.path("outcomes", "secondary", "7-2-anticoagulation-concurrent-intermediate-summary"))
make_summary_table(acs_itt_concurc2_dat)
Anticoagulation intervention Patients Known Deaths Hospitalised WHO, Median (Q1, Q3)
Low-dose 610 596 19 (3%) 7 (1%) 1 (1, 2)
Intermediate-dose 613 603 15 (2%) 5 (1%) 1 (1, 2)
Overall 1223 1199 34 (3%) 12 (1%) 1 (1, 2)
Sample distribution
pdat <- acs_itt_concurc2_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2"),
      labels = c("Low", "Intermediate")),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = who), colour = NA) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-2-descriptive-concurrent-intermediate.pdf"), p, height = 2.5, width = 4)
p

Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc2_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc2_nona_dat$D28_who)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 75136,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(2) %**% mdrws$beta[1])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt), 
              setNames(exp(mdrws$beta[2:4]), c("Age \u2265 60", "Australia/New Zealand", "Nepal")))
Posterior contrast
p <- plot_or_densities(mdrws$compare)
p

Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-2-primary-model-acs-itt-concurrent-intermediate-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.83 (0.64, 1.08) 0.84 (0.11) 0.92
Age ≥ 60 2.77 (2.11, 3.65) 2.80 (0.39) 0.00
Australia/New Zealand 2.21 (1.51, 3.23) 2.25 (0.44) 0.00
Nepal 0.56 (0.31, 0.97) 0.59 (0.17) 0.98
Posterior summaries for model parameters.
Code
mfit$summary()
# A tibble: 24 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -8.80e+2 -8.80e+2 2.43e+0 2.30e+0 -8.85e+2 -8.77e+2  1.00    8198.
 2 beta_raw[… -1.29e-1 -1.29e-1 9.33e-2 9.34e-2 -2.83e-1  2.32e-2  1.00   32680.
 3 beta_raw[…  4.08e-1  4.08e-1 5.53e-2 5.45e-2  3.18e-1  4.99e-1  1.00   21528.
 4 beta_raw[…  7.94e-1  7.95e-1 1.95e-1 1.94e-1  4.73e-1  1.11e+0  1.00   27136.
 5 beta_raw[… -5.77e-1 -5.71e-1 2.88e-1 2.86e-1 -1.06e+0 -1.12e-1  1.00   32058.
 6 p[1]        7.91e-1  7.91e-1 1.45e-2 1.45e-2  7.66e-1  8.14e-1  1.00   19596.
 7 p[2]        1.84e-1  1.84e-1 1.30e-2 1.30e-2  1.63e-1  2.06e-1  1.00   21477.
 8 p[3]        1.22e-3  1.05e-3 8.10e-4 7.06e-4  2.47e-4  2.76e-3  1.00   27087.
 9 p[4]        2.29e-3  2.11e-3 1.13e-3 1.07e-3  8.16e-4  4.42e-3  1.00   27665.
10 p[5]        1.20e-3  1.03e-3 8.11e-4 7.02e-4  2.43e-4  2.78e-3  1.00   29866.
# … with 14 more rows, and 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9741940 1.0757049 1.0154633 1.0669753 0.9840496 1.0026432 0.9969708
[8] 0.9689346
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])

Low-dose with aspirin

  • Set: ACS-ITT-aspirin
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc3_dat, "latex"), 
  file.path("outcomes", "secondary", "7-2-anticoagulation-concurrent-stdaspirin-summary"))
make_summary_table(acs_itt_concurc3_dat)
Anticoagulation intervention Patients Known Deaths Hospitalised WHO, Median (Q1, Q3)
Low-dose 299 291 11 (4%) 5 (2%) 1 (1, 2)
Intermediate-dose 298 293 12 (4%) 3 (1%) 1 (1, 2)
Low-dose with aspirin 283 281 10 (4%) 5 (2%) 1 (1, 2)
Overall 880 865 33 (4%) 13 (2%) 1 (1, 2)
Sample distribution
pdat <- acs_itt_concurc3_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3"),
      labels = c("Low", "Intermediate", "Low\nwith aspirin")),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = who)) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-2-descriptive-concurrent-stdaspirin.pdf"), p, height = 2.5, width = 6)
p

Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc3_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc3_nona_dat$D28_who)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1)
)
snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 135356,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Low with aspirin" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:4]), c("Age \u2265 60", "Australia/New Zealand")))
Posterior contrast
plot_or_densities(mdrws$compare)

Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-2-primary-model-acs-itt-concurrent-stdaspirin-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.86 (0.61, 1.22) 0.87 (0.16) 0.80
Low with aspirin 0.73 (0.51, 1.05) 0.74 (0.14) 0.96
Age ≥ 60 2.10 (1.53, 2.88) 2.13 (0.34) 0.00
Australia/New Zealand 0.95 (0.52, 1.68) 0.99 (0.30) 0.56
Posterior summaries for model parameters.
Code
mfit$summary()
# A tibble: 24 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -7.14e+2 -7.13e+2 2.45    2.30    -7.18e+2 -7.10e+2  1.00    7984.
 2 beta_raw[… -1.17e-1 -1.16e-1 0.131   0.132   -3.33e-1  9.73e-2  1.00   37321.
 3 beta_raw[…  1.91e-1  1.91e-1 0.126   0.124   -1.81e-2  3.98e-1  1.00   39109.
 4 beta_raw[…  2.97e-1  2.97e-1 0.0643  0.0639   1.91e-1  4.03e-1  1.00   24622.
 5 beta_raw[… -5.40e-2 -4.95e-2 0.300   0.297   -5.59e-1  4.26e-1  1.00   34330.
 6 p[1]        7.43e-1  7.43e-1 0.0178  0.0177   7.14e-1  7.72e-1  1.00   23448.
 7 p[2]        2.14e-1  2.14e-1 0.0156  0.0154   1.89e-1  2.40e-1  1.00   26429.
 8 p[3]        2.08e-3  1.78e-3 0.00138 0.00121  4.26e-4  4.79e-3  1.00   31594.
 9 p[4]        3.89e-3  3.58e-3 0.00189 0.00175  1.39e-3  7.44e-3  1.00   33907.
10 p[5]        2.97e-3  2.66e-3 0.00167 0.00150  8.59e-4  6.17e-3  1.00   31434.
# … with 14 more rows, and 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9850497 1.0408767 0.9642583 1.0151688 1.0564695 0.9240525 0.9833757
[8] 0.9764285
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])

Therapeutic dose

  • Set: ACS-ITT-therapeutic
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc4_dat, "latex"), 
  file.path("outcomes", "secondary", "7-2-anticoagulation-concurrent-therapeutic-summary"))
make_summary_table(acs_itt_concurc4_dat)
Anticoagulation intervention Patients Known Deaths Hospitalised WHO, Median (Q1, Q3)
Low-dose 79 75 3 (4%) 1 (1%) 1 (1, 2)
Intermediate-dose 65 63 1 (2%) 1 (2%) 1 (1, 2)
Therapeutic-dose 50 50 6 (12%) 1 (2%) 1 (1, 2)
Overall 194 188 10 (5%) 3 (2%) 1 (1, 2)
Sample distribution
pdat <- acs_itt_concurc4_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C4"),
      labels = c("Low", "Intermediate", "Therapeutic")),
    who = ordered(D28_who, levels = 1:8), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = who)) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("WHO scale", option = "A", direction = -1) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-2-descriptive-concurrent-therapeutic.pdf"), p, height = 2.5, width = 6)
p

Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc4_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc4_nona_dat$D28_who)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1, 1)
)

snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 49135,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Therapeutic vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Therapeutic" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:5]), c("Age \u2265 60", "India", "Australia/New Zealand")))
Posterior contrasts
plot_or_densities(mdrws$compare)

Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/secondary/7-2-primary-model-acs-itt-concurrent-therapeutic-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.49 (0.22, 1.03) 0.53 (0.21) 0.97
Therapeutic 1.08 (0.50, 2.32) 1.16 (0.47) 0.42
Age ≥ 60 3.27 (1.68, 6.45) 3.48 (1.24) 0.00
India 0.99 (0.37, 2.57) 1.11 (0.58) 0.51
Australia/New Zealand 3.87 (2.02, 7.63) 4.12 (1.45) 0.00
Posterior summaries for model parameters.
Code
mfit$summary()
# A tibble: 22 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -1.67e+2 -1.67e+2 2.38    2.23    -1.71e+2 -1.64e+2  1.00    8588.
 2 beta_raw[…  5.59e-1  5.59e-1 0.292   0.290    8.38e-2  1.04e+0  1.00   22664.
 3 beta_raw[…  2.61e-1  2.61e-1 0.270   0.268   -1.84e-1  7.08e-1  1.00   20274.
 4 beta_raw[…  4.75e-1  4.74e-1 0.137   0.136    2.52e-1  7.01e-1  1.00   14420.
 5 beta_raw[… -1.40e-2 -9.69e-3 0.496   0.491   -8.38e-1  7.88e-1  1.00   20324.
 6 beta_raw[…  1.36e+0  1.35e+0 0.338   0.339    8.08e-1  1.92e+0  1.00   14624.
 7 p[1]        8.52e-1  8.55e-1 0.0367  0.0359   7.86e-1  9.07e-1  1.00   11710.
 8 p[2]        1.25e-1  1.23e-1 0.0306  0.0301   7.91e-2  1.80e-1  1.00   12950.
 9 p[3]        2.39e-3  1.72e-3 0.00228 0.00161  2.02e-4  6.85e-3  1.00   18472.
10 p[4]        2.25e-3  1.63e-3 0.00211 0.00151  1.86e-4  6.42e-3  1.00   15179.
# … with 12 more rows, and 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9336335 0.9877803 1.0642128 0.9765569 0.9757024 1.0348876 1.0237602
[8] 1.0054517
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])

Sensitivity: Treatment contrast prior

Fit primary model
res <- fit_primary_model(acs_itt_nona_dat, contr.treatment, save = FALSE)
Odds ratio summary table
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.79 (0.60, 1.05) 0.80 (0.12) 0.94
Low with aspirin 0.75 (0.52, 1.06) 0.76 (0.14) 0.95
Therapeutic 1.45 (0.73, 2.84) 1.54 (0.55) 0.14
Ineligible aspirin 1.71 (0.77, 3.70) 1.85 (0.77) 0.09
Age ≥ 60 2.48 (1.90, 3.25) 2.50 (0.35) 0.00
Australia/New Zealand 1.41 (0.42, 4.32) 1.66 (1.05) 0.28
Nepal 0.75 (0.23, 2.46) 0.90 (0.64) 0.70
Posterior summaries for model parameters (fixed-effects).
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch),
  exp(res$drws$gamma_site),
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
p

Code
p <- plot_or_densities(res$drws$compare)
p